library(haven) # 读取NHANES XPT文件的包
library(plyr) # 用于数据处理
library(dplyr) # 用于数据处理
library(arsenal) # 用于数据快速预览
library(survey) # 用于加权情况下的分析

setwd("D:\\NHANES DATA") # 设置工作目录

demo_temp = read_xpt("1673147009143.xpt") # DEMO-人口学数据提取：2017-2018
View(demo_temp)

# list.files()
demo_j = read_xpt("2017-2018\\Demographics\\demo_j.xpt") # DEMO-人口学数据提取：2017-2018
# colnames(demo_j)

# dim(demo_j) # [1] 9254   46
# colnames(demo_g) # 输出列名，即变量名称
# View(demo_j)
# table_1 = tableby( ~ RIDAGEYR + factor(RIAGENDR) + RIDRETH3 + DMDEDUC2 + INDFMPIR, data = demo_g)
# summary(table_1 , text = TRUE)
demo_data = demo_j[ , c("SEQN" , "RIDAGEYR" , "RIAGENDR" , "RIDRETH3" , "DMDEDUC2")] # 只有其中的几个变量
# View(demo_data)

paq_j = read_xpt("2017-2018\\Questionnaire\\paq_j.xpt") # 1.2 PAQ-运动数据提取
paq_j_data = paq_j[ , c("SEQN" , "PAQ650" , "PAQ665")]

bmx_j = read_xpt("2017-2018\\Examination\\bmx_j.xpt") # 1.3 BMI 数据提取
bmx_j_data = bmx_j[ , c("SEQN" , "BMXBMI")]

smq_j = read_xpt("2017-2018\\Questionnaire\\smq_j.xpt") # 1.4 SMQ-吸烟数据提取
smq_j_data = smq_j[ , c("SEQN" , "SMQ020" , "SMQ040" , "SMQ050Q")]

# 1.5 糖尿病诊断
ghb_j = read_xpt("2017-2018\\Laboratory\\ghb_j.xpt") # 1.5.1 糖化血红蛋白A1C（%）
ghb_j_data = ghb_j[ , c("SEQN" , "LBXGH")]

glu_j = read_xpt("2017-2018\\Laboratory\\glu_j.xpt") # 1.5.2 空腹血糖（mg/dl）
glu_j_data = glu_j[ , c("SEQN" , "LBXGLU")]

diq_j = read_xpt("2017-2018\\Questionnaire\\diq_j.xpt") # 1.5.3 是否有医生告知您患有糖尿病、现在是否使用胰岛素
diq_j_data = diq_j[ , c("SEQN" , "DIQ010" , "DIQ050")]

diabetes = plyr::join_all(list(ghb_j_data, glu_j_data, diq_j_data))
# dim(diabetes) # [1] 6401    5
diabetes$diabetes_index = ifelse(diabetes$LBXGH >= 6.5 | diabetes$LBXGLU >= 126 | diabetes$DIQ010 ==1 | diabetes$DIQ050 == 1 , 1 , NA)
# table(diabetes_index) # 1112 
diabetes_data = diabetes[ , c("SEQN" , "diabetes_index")]

alt_j = read_xpt("2017-2018\\Laboratory\\biopro_j.xpt") # 1.6 丙氨酸转氨酶
alt_j_data = biopro_j[ , c("SEQN" , "LBXSATSI")]

vic_j = read_xpt("2017-2018\\Laboratory\\vic_j.xpt") # 1.7 血清维生素C的含量
vic_j_data = vic_j[ , c("SEQN" , "LBDVICSI" , "LBXVIC")]

hdl_j = read_xpt("2017-2018\\Laboratory\\hdl_j.xpt") # 1.8 高密度脂蛋白
hdl_j_data = hdl_j[ , c("SEQN" , "LBDHDDSI" , "LBDHDD")]

lux_j = read_xpt("2017-2018\\Examination\\lux_j.xpt") # 1.9 弹性成像检查
lux_j_data = lux_j[ , c("SEQN" , "LUAXSTAT" , "LUXSMED" , "LUXCAPM")]
# View(lux_j_data)

dr1tot_j = read_xpt("2017-2018\\Dietary\\dr1tot_j.xpt")
dr2tot_j = read_xpt("2017-2018\\Dietary\\dr2tot_j.xpt")
dr_data_file = merge(dr1tot_j , dr2tot_j) # merge这个是拼接，默认使用两个数据框中相同列名称进行拼接。
alco_data = dr_data_file[ , c("SEQN" , "DR1TALCO" , "DR2TALCO")] # 只有其中的几个变量
# dim(alco_data) # [1] 8704    3

# 1.11 其他原因导致的肝脏疾病
hepbd_j = read_xpt("2017-2018\\Laboratory\\hepbd_j.xpt") # 乙型肝炎表面抗原
hepbd_j_data = hepbd_j[ , c("SEQN" , "LBDHBG")]

hepc_j = read_xpt("2017-2018\\Laboratory\\hepc_j.xpt") # 丙肝抗体\丙肝RNA
hepc_j_data = hepc_j[ , c("SEQN" , "LBDHCI" , "LBXHCR")]

mcq_j = read_xpt("2017-2018\\Questionnaire\\mcq_j.xpt") # 肝癌、自身免疫性肝炎
mcq_j_data = mcq_j[ , c("SEQN" , "MCQ510E" , "MCQ230A" , "MCQ230B" , "MCQ230C")]

weight_data = demo_j[ , c("SEQN" , "WTMEC2YR")] # 权重的计算方式，需要按 Cycle 进行平均，这个仅有1个 Cycle，无需做处理

survey_design_data = demo_j[ , c("SEQN" , "SDMVPSU" , "SDMVSTRA")] # 2.2 复杂抽样的其他变量-DEMO

output <- plyr::join_all(list(demo_data , paq_j_data , bmx_j_data , smq_j_data , ghb_j_data , glu_j_data , diq_j_data , alt_j_data , vic_j_data , hdl_j_data , 
                              lux_j_data , alco_data , hepbd_j_data , hepc_j_data , mcq_j_data , weight_data , survey_design_data , diabetes_data), by="SEQN", type="left")
# View(output)
 dim(output) # [1] 9254   34

# length(which(is.na(output$LBDVICSI))) # [1] 2514 #  4.1缺少血清维生素 C 数据(n = 2514)
data_vc_exist = output[ - which(is.na(output$LBDVICSI)) , ] # 已减去缺少血清维生素 C 的数据
# dim(data_vc_exist) # [1] 6740   34
# colnames(data_vc_exist)

# 4.2 排除 LUX 数据缺失以及不是 NAFLD 的人员 弹性成像检查状态-LUAXSTAT; 可控衰减参数中位数-LUXCAPM; 中位硬度-LUXSMED
cap_exclude = which(data_vc_exist$LUAXSTAT != 1 | is.na(data_vc_exist$LUXCAPM) | is.na(data_vc_exist$LUXSMED) | data_vc_exist$LUXCAPM < 248)
# length(cap_exclude) # [1] 3945
data_vc_exist_cap_exist = data_vc_exist[ - cap_exclude , ]
# dim(data_vc_exist_cap_exist) # [1] 2795   34

# 4.3 排除其他原因导致的肝脏疾病 & 重要协变量的缺失
liver_index = which(data_vc_exist_cap_exist$MCQ230A == 22 | data_vc_exist_cap_exist$MCQ230B == 22 | data_vc_exist_cap_exist$MCQ230C == 22) # 肝癌 MCQ230A-C
# length(liver_index) # [1] 1
autoimm_hepa_index = which(data_vc_exist_cap_exist$MCQ510E == 5)
# length(autoimm_hepa_index) # [1] 4
hepc_index = which(data_vc_exist_cap_exist$LBDHCI == 1 | data_vc_exist_cap_exist$LBXHCR == 1)
# length(hepc_index) # [1] 40

hepbd_index = which(data_vc_exist_cap_exist$LBDHBG == 1)
# length(hepbd_index) # [1] 15


# dim(data_vc_exist_cap_exist)
# View(data_vc_exist_cap_exist[ , c("SEQN" , "DR1TALCO" , "DR2TALCO")])
# 如果两天均有值，则使用两天的平均，如果第二天为空，则使用第一天的值
data_vc_exist_cap_exist$total_dr = ifelse(is.na(data_vc_exist_cap_exist$DR2TALCO) , data_vc_exist_cap_exist$DR1TALCO , (data_vc_exist_cap_exist$DR1TALCO + data_vc_exist_cap_exist$DR2TALCO) / 2) 
# View(data_vc_exist_cap_exist[ , c("SEQN" , "total_dr" , "DR1TALCO" , "DR2TALCO")])
# 是否过量饮酒，1为过量，0为不过量，男大于20，女大于10
data_vc_exist_cap_exist$total_drink = ifelse(data_vc_exist_cap_exist$RIAGENDR == 1 , ifelse(data_vc_exist_cap_exist$total_dr >= 20 , 1 , 0) , ifelse(data_vc_exist_cap_exist$total_dr >= 10 , 1 , 0)) 
# View(data_vc_exist_cap_exist[ , c("SEQN" , "total_drink" , "RIAGENDR" , "total_dr")])
# dim(data_vc_exist_cap_exist)

excessive_alco_index = which(data_vc_exist_cap_exist$total_drink == 1) # 过度饮酒
# length(excessive_alco_index) # [1] 371

other_index = c(autoimm_hepa_index , hepbd_index , hepc_index , liver_index , excessive_alco_index)

data_vc_exist_cap_exist_not_other = data_vc_exist_cap_exist[ - other_index , ]
# dim(data_vc_exist_cap_exist_not_other) # [1] 2374   36

# covariate 上的缺失数据  体重指数和是否吸烟至少100支这两项，如果添加上，数据为1917条，与论文中的不相符，所以先给注掉。
paper_data = subset.data.frame(data_vc_exist_cap_exist_not_other , 
                    (!is.na(LBDVICSI)) & # 维C
                    (!is.na(LBXSATSI)) & # 丙氨酸转氨酶
                    (!is.na(LBDHDDSI)) & # 高密度脂蛋白
                   # (!is.na(BMXBMI)) & # 体重指数
                   # (!is.na(SMQ020)) & # 是否吸烟至少100支
                    (!is.na(RIDRETH3)) & # 种族
                    (!is.na(LUXCAPM)) & # 可控衰减参数中位数
                    (!is.na(LBXGH)) & # 糖化血红蛋白A1C(%)
                    (!is.na(total_dr)) & # 两天平均的饮酒量
                    (!is.na(RIDAGEYR)) & # 年龄
                    (!is.na(RIAGENDR)) & # 性别
                    (!is.na(RIDRETH3)) & # 种族
                    (!is.na(DMDEDUC2))  # 教育程度
)
# dim(paper_data) # [1] 1926   36

# 5. 变量结果 Summary
nhanes_design <- svydesign(data = paper_data[which(paper_data$WTMEC2YR > 0) , ] , id = ~ SDMVPSU , strata = ~ SDMVSTRA , weights = ~ WTMEC2YR , nest = TRUE)
summary(nhanes_design, text=TRUE)
# View(paper_data[ , c("WTMEC2YR" , "SDMVPSU" , "SDMVSTRA")])

# 2.1 变量（含加权）Sample
svymean( ~ RIDAGEYR + factor(RIAGENDR) , nhanes_design, na.rm=TRUE) # 和 survey total 是不一致的

# 2.2 分类变量-count结果 Sample
tab1 = tableby( ~ factor(RIAGENDR) + factor(diabetes_index) , data = paper_data) 
summary(tab1 , text = TRUE)

#### 6. 课后作业 ####
##### 6.1  衍生 diabetes 指标#####
# 不参考原始代码，重新衍生 Diabetes，目前我们只将有 diabetes 标记为1，没有的则是 NA
# 将它衍生为 Yes/No 的取值
#----------------------------------- 从这里开始写代码
# ----糖尿病
paper_data$diabetes_data = ifelse(is.na(paper_data$diabetes_index) | paper_data$diabetes_index == 0 , "No" , "Yes")
# View(paper_data[ , c("diabetes_index" , "SEQN" , "diabetes_data")]) # diabetes_index这里的只有1和NA，diabetes_data这里的给改成了Yes或No
tab_diabetes = tableby( ~ factor(diabetes_data)  , data = paper_data) 
summary(tab_diabetes) # No：1345 (69.8%)；Yes：581 (30.2%) 与论文差距较大

##### 6.2 衍生其他分类变量 #####
# 衍生 Education level、BMI group、 Recreational physical activity、Smoking status 的分类变量，
# 命名为 edu.group, bmi.group, recre.phy.group, smoke.group
#----------------------------------- 从这里开始写代码
# ----教育程度
# View(paper_data[ , c("SEQN" , "DMDEDUC2")])
paper_data$edu_group = ifelse(paper_data$DMDEDUC2 <= 3 , "<=High school" , ">High school")
# View(paper_data[ , c("DMDEDUC2" , "SEQN" , "edu_group")])
tab_edu = tableby( ~ factor(edu_group)  , data = paper_data) 
summary(tab_edu) # <=High school：861 (44.7%)；>High school：1065 (55.3%) 与论文百分比差距较大

# ----体重指数分组
# View(paper_data[ , c("SEQN" , "BMXBMI")])
b = paper_data$BMXBMI
paper_data$BMI_group = ifelse((is.na(b) | b < 25), "<25" , ifelse(b >= 25 & b <= 30 , "25-30", ifelse(b > 30 , ">30" , "NA")))
# View(paper_data[ , c("BMXBMI" , "SEQN" , "BMI_group")])
tab_bmi = tableby( ~ factor(BMI_group)  , data = paper_data) 
summary(tab_bmi) # <25：215 (11.2%)；25-30：611 (31.7%)；>30：1100 (57.1%) 与论文差距较大

# ----娱乐体育活动
# View(paper_data[ , c("SEQN" , "PAQ665")])
paper_data$recre_phy_group = ifelse(paper_data$PAQ665 == 1 , "Yes" , "No")
# View(paper_data[ , c("PAQ665" , "SEQN" , "recre_phy_group")])
tab_recre_phy = tableby( ~ factor(recre_phy_group)  , data = paper_data) 
summary(tab_recre_phy) # No：1208 (62.7%)；Yes：718 (37.3%)

# ----吸烟状态
# View(paper_data[ , c("SEQN" , "SMQ020" , "SMQ040")])
# SMQ020 - 一生中至少吸过 100 支香烟 2：不是（按从不吸烟Never）1：是的（SMQ040 - 你现在抽烟吗  3：从不，回答3按以前的吸烟者Former，回答其它按当前的吸烟者Current）
paper_data$smoke_group = ifelse((paper_data$SMQ020 == 2) , "Never", ifelse((paper_data$SMQ040 == 3) , "Former" , "Current"))
# View(paper_data[ , c("SMQ020" , "SMQ040" , "smoke_group")])
tab_smoke = tableby( ~ factor(smoke_group)  , data = paper_data) 
summary(tab_smoke) # Current：272 (14.1%)；Former：533 (27.7%) ；Never：1121 (58.2%) 与论文差距较大

##### 6.3 复现 Table1-1 #####
# 复现 Table1 第一列离散变量的结果，使用 tableby 函数输出分类变量在全部人群中的指标(count)
#----------------------------------- 从这里开始写代码 （说明：教育程度、体重指数、娱乐性体育活动、吸烟状态上面已有，下面就不重复写了）

# ----年龄的平均值和标准差
mean(paper_data$RIDAGEYR , na.rm = TRUE) # [1] 54.06023
sd(paper_data$RIDAGEYR, na.rm = TRUE) # [1] 16.02085

# ----性别
# View(paper_data[ , c("SEQN" , "RIAGENDR")])
paper_data$gender = ifelse(paper_data$RIAGENDR == 1 , "Male" , "Female")
# View(paper_data[ , c("SEQN" , "RIAGENDR" , "gender")])
tab_gender = tableby( ~ factor(gender)  , data = paper_data) 
summary(tab_gender) # Female：938 (48.7%)；Male：988 (51.3%)

# ----种族
# View(paper_data[ , c("SEQN" , "RIDRETH3")])
r = paper_data$RIDRETH3
paper_data$race = ifelse((r == 1), "Mexican American" , ifelse(r == 2 , "Other Hispanic",  ifelse(r == 3 , "Non-Hispanic white",  ifelse(r == 4 , "Non-Hispanic black", "Other race"))))
# View(paper_data[ , c("SEQN" , "RIDRETH3" , "race")])
tab_race = tableby( ~ factor(race)  , data = paper_data) 
summary(tab_race) # Mexican American：340 (17.7%)；Non-Hispanic black：362 (18.8%)；Non-Hispanic white：673 (34.9%)；Other Hispanic：191 (9.9%)；Other race：360 (18.7%)

# ----血清维生素C（注意，这里使用的是LBXVIC，因为这个的单位是mg/dL，与论文中的相同）
# View(paper_data[ , c("SEQN" , "LBXVIC")])
mean(paper_data$LBXVIC , na.rm = TRUE) # [1] 0.8532336
sd(paper_data$LBXVIC, na.rm = TRUE) # [1] 0.4748143

# ----丙氨酸转氨酶
# View(paper_data[ , c("SEQN" , "LBXSATSI")])
mean(paper_data$LBXSATSI , na.rm = TRUE) # [1] 24.1298
sd(paper_data$LBXSATSI, na.rm = TRUE) # [1] 15.92155

# ----高密度脂蛋白（注意，这里使用的是LBDHDD，因为这个的单位是mg/dL，与论文中的相同）
# View(paper_data[ , c("SEQN" , "LBDHDD")])
mean(paper_data$LBDHDD , na.rm = TRUE) # [1] 48.5026
sd(paper_data$LBDHDD, na.rm = TRUE) # [1] 13.16144

# ----中位可控衰减参数
# View(paper_data[ , c("SEQN" , "LUXCAPM")])
mean(paper_data$LUXCAPM , na.rm = TRUE) # [1] 307.163
sd(paper_data$LUXCAPM, na.rm = TRUE) # [1] 40.58216

##### 6.4 复现 Table1-2 #####
# 复现 Table1 第一列连续变量/离散变量的结果，使用 svymean 函数输出分类变量在全部人群中的指标
#----------------------------------- 从这里开始写代码 
svymean( ~ paper_data$LBXVIC , nhanes_design, na.rm=TRUE) # 和 survey total 是不一致的
#.................................这个不知如何作..............................................



##### 6.5 复现 Table1-3 #####
# 衍生 significant fibrosis 指标，将人群分为 significant fibrosis == yes & no
# 命名为 fibrosis.group
#----------------------------------- 从这里开始写代码 
# ----中值刚度
# View(paper_data[ , c("SEQN" , "LUXSMED")])
paper_data$fibrosis_group = ifelse(paper_data$LUXSMED >= 8.2 , "Yes" , "No")
# View(paper_data[ , c("SEQN" , "LUXSMED" , "fibrosis_group")])
tab_fibrosis_group = tableby( ~ factor(fibrosis_group)  , data = paper_data) 
summary(tab_fibrosis_group) # No：1661 (86.2%)；Yes：265 (13.8%)


##### 6.6 复现 Table1-4 #####
# 复现 Table1 第2列 & 第3列根据不同 fibrosis.group 的统计结果（离散变量-count数）
# 使用 tableby 函数输出分类变量在significant fibrosis 不同分层中的结果
#----------------------------------- 从这里开始写代码 

fibrosis_yes = paper_data[which(paper_data$fibrosis_group == "Yes") , ]
fibrosis_no = paper_data[which(paper_data$fibrosis_group == "No") , ]

# dim(fibrosis_yes) # [1] 265  47
# dim(fibrosis_no) # [1] 1661   47
# dim(paper_data) # [1] 1926   47

# ~~~~年龄的平均值和标准差
mean(fibrosis_yes$RIDAGEYR , na.rm = TRUE) # [1] 56.28302
sd(fibrosis_yes$RIDAGEYR, na.rm = TRUE) # [1] 14.6002
mean(fibrosis_no$RIDAGEYR , na.rm = TRUE) # [1] 53.7056
sd(fibrosis_no$RIDAGEYR, na.rm = TRUE) # [1] 16.21191

# ~~~~性别
summary(tableby( ~ factor(gender)  , data = fibrosis_yes)) # Female：115 (43.4%)；Male：150 (56.6%)
summary(tableby( ~ factor(gender)  , data = fibrosis_no)) # Female：823 (49.5%)；Male：838 (50.5%)

# ~~~~种族
summary(tableby( ~ factor(race)  , data = fibrosis_yes)) # Mexican American：50 (18.9%)；Non-Hispanic black：48 (18.1%)；Non-Hispanic white：97 (36.6%)；Other Hispanic：30 (11.3%)；Other race：40 (15.1%)
summary(tableby( ~ factor(race)  , data = fibrosis_no)) # Mexican American：290 (17.5%)；Non-Hispanic black：314 (18.9%)；Non-Hispanic white：576 (34.7%)；Other Hispanic：161 (9.7%)；Other race：320 (19.3%)

# ~~~~教育程度
summary(tableby( ~ factor(edu_group)  , data = fibrosis_yes)) # <=High school：128 (48.3%)；>High school：137 (51.7%)
summary(tableby( ~ factor(edu_group)  , data = fibrosis_no)) # <=High school：733 (44.1%)；>High school：928 (55.9%)

# ~~~~体重指数分组
summary(tableby( ~ factor(BMI_group)  , data = fibrosis_yes)) # <25：9 (3.4%)；25-30：37 (14.0%)；>30：219 (82.6%)
summary(tableby( ~ factor(BMI_group)  , data = fibrosis_no)) # <25：206 (12.4%)；25-30：574 (34.6%)；>30：881 (53.0%)

# ~~~~娱乐性体育活动 
summary(tableby( ~ factor(recre_phy_group)  , data = fibrosis_yes)) # No：178 (67.2%)；Yes：87 (32.8%)
summary(tableby( ~ factor(recre_phy_group)  , data = fibrosis_no)) # No：1030 (62.0%)；Yes：631 (38.0%)

# ~~~~吸烟状态
summary(tableby( ~ factor(smoke_group)  , data = fibrosis_yes)) # Current：31 (11.7%)；Former：92 (34.7%) ；Never：142 (53.6%)
summary(tableby( ~ factor(smoke_group)  , data = fibrosis_no)) # Current：241 (14.5%)；Former：441 (26.6%) ；Never：979 (58.9%)

# ~~~~糖尿病
summary(tableby( ~ factor(diabetes_data)  , data = fibrosis_yes)) # No：126 (47.5%)；Yes：139 (52.5%)
summary(tableby( ~ factor(diabetes_data)  , data = fibrosis_no)) # No：1219 (73.4%)；Yes：442 (26.6%)

# ~~~~血清维生素C（注意，这里使用的是LBXVIC，因为这个的单位是mg/dL，与论文中的相同）
mean(fibrosis_yes$LBXVIC , na.rm = TRUE) # [1] 0.7388038
sd(fibrosis_yes$LBXVIC, na.rm = TRUE) # [1] 0.4849571
mean(fibrosis_no$LBXVIC , na.rm = TRUE) # [1] 0.8714901
sd(fibrosis_no$LBXVIC, na.rm = TRUE) # [1] 0.4707573

# ~~~~丙氨酸转氨酶
mean(fibrosis_yes$LBXSATSI , na.rm = TRUE) # [1] 30.4566
sd(fibrosis_yes$LBXSATSI, na.rm = TRUE) # [1] 23.40533
mean(fibrosis_no$LBXSATSI , na.rm = TRUE) # [1] 23.12041
sd(fibrosis_no$LBXSATSI, na.rm = TRUE) # [1] 14.12206

# ~~~~高密度脂蛋白（注意，这里使用的是LBDHDD，因为这个的单位是mg/dL，与论文中的相同）
mean(fibrosis_yes$LBDHDD , na.rm = TRUE) # [1] 44.82642
sd(fibrosis_yes$LBDHDD, na.rm = TRUE) # [1] 12.6566
mean(fibrosis_no$LBDHDD , na.rm = TRUE) # [1] 49.0891
sd(fibrosis_no$LBDHDD, na.rm = TRUE) # [1] 13.14912

# ~~~~中位可控衰减参数
mean(fibrosis_yes$LUXCAPM , na.rm = TRUE) # [1] 334.7396
sd(fibrosis_yes$LUXCAPM, na.rm = TRUE) # [1] 43.19209
mean(fibrosis_no$LUXCAPM , na.rm = TRUE) # [1] 302.7634
sd(fibrosis_no$LUXCAPM, na.rm = TRUE) # [1] 38.37138

#####  6.7 复现 Table1-5  #####
# 复现 Table1 第2列 & 第3列根据不同 fibrosis.group 的统计结果（连续变量/离散变量）
# 使用 svymean 函数输出分类变量在全部人群中的指标
#----------------------------------- 从这里开始写代码 
#.................................这个不知如何作..............................................











